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ABSTRACT 

We present results from a large volume simulation of Hydrogen reionization. We combine 3d radia- 
tive transfer calculations and an N-body simulation, describing structure formation in the intergalactic 
medium (IGM), to detail the growth of HII regions around high redshift galaxies. Our N-body sim- 
ulation tracks 1024 3 dark matter particles, in a cubical box of co-moving side length L^ ox = 65.6 
Mpc h . This large volume allows us to accurately characterize the size distribution of HII regions 
throughout most of the reionization process. At the same time, our simulation resolves many of the 
small galaxies likely responsible for reionization. It confirms a picture anticipated by analytic models: 
HII regions grow collectively around highly-clustered sources, and have a well-defined characteristic 
size, which evolves from a sub-Mpc scale at the beginning of reionization to R > 10 co-moving Mpc 
towards the end. We show that in order to obtain this qualitative picture, source resolution must 
not be sacrificed at too great a level. We present a detailed statistical description of our results, and 
compare them with a numerical hybrid scheme based on the analytic model by Furlanetto, Zaldar- 
riaga, and Hernquist. This model associates HII regions with large-scale overdensities and is based 
on the excursion set formalism. We find that the analytic calculation reproduces the size distribution 
of HII regions, the power spectrum of the ionization field, and the 21 cm power spectrum of the 
full radiative transfer simulation remarkably well. The ionization field from the radiative transfer 
simulation, however, has more small scale structure than the analytic calculation, owing to Poisson 
scatter in the simulated abundance of galaxies on small scales. We propose and validate a simple 
scheme to incorporate this scatter into our calculations. Our results suggest that analytic calculations 
are sufficiently accurate to aid in predicting and interpreting the results of future 21 cm surveys. In 
particular, our fast numerical scheme is useful for forecasting constraints from future 21 cm surveys, 
and in constructing mock surveys to test data analysis procedures. 

Subject headings: cosmology: theory - intergalactic medium - large scale structure of universe 



1. INTRODUCTION 

The epoch of reionization (EOR) is a pivotal stage in 
the process of cosmological structure formation, marking 
the birth of the first luminous objects, a key landmark as 
the universe transforms from the relatively smooth state 
probed by the cosmic microwave background (CMB), to 
its present day complexity. Our current observational 
constraints on reionization come from Lya forest ab- 
sorption spectra towards high redshift quasars and con- 
straints on the ev olution of the ionizing background (e.g. 
Fan et al. 2005), from measurements of the high red- 
shift galaxy lumin osity function from narrow- band Lya- 



meager. Quasar absorption spectra are limited in part 
by the high Lya absorption cross section: by z ~ 6, 
even a highly ionized IGM completely absorbs quasar 
flux in the Lya forest. The constraints from narrow- 



emission searches (Malhotra & Rhoads 
measurements of trie large scale (JMB 



tion (Page et al.|2006 Spergel et al 



2006) 



2005], and from 
mode polariza- 
The claimed 



size of Hll regions surrounding individual quasars has 
also been used to infer limits on the _neutral fraction 

While valu- 

., i ~as its limita- ; 

tions, and some of the current constraints are relatively (e.g. IMaaau et al.||HJiJo 



( |Mesinger et al |2004| [Wyithe et aL|2004| ). 
able, each of these observational probes n. 



band Lya search es are subtle to interpret (e.g. Furlan- 
etto et al.|2006b I , and restricted to narrow redshift win- 
dows around z = 5.7 and z = 6.5, where Lya falls in the 
observed optical ba nd, and avoids cont amination from 
bright sky lines (e.g. |Rhoads et al.|2003[ ). The CMB po- 
larization measurements constrain only an integral over 
the ionization history, an d are potentially sensitive to 
foreg round contamination ( Kogut et al.|2003 Page et al 
20061. 

The study of reionization may be revolutionized by 
future experiments aimed at detecting 21 cm emission 
from the high redshift IGM. These experiments should 
provide three-dimensional information regarding the dis- 
tribution of high redshift neutral hydrogen, constraining 
the t opology of reionization, and its redshift evolution 
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Zaldarriaga et al. 2004 1. Sev- 
eral low frequency radio telescopes are presently ramp- 
ing up to dete ct this signal: the M ileura Wide Field 
Array (MWA) ( |Bowman et al.||2006[> R the P rimeavAl 
Structure Telescope (PAST) (|Pen et al.||2004[), and the 



http://web.haystack.mit.edu/arrays/MWA/ 



Low Frequency Array (LOFAR) f\ while another sec- 
ond generation experiment, the Square Kilometer Array 
(SKArl is in the planning stage. These measurements 
will be dominated by foreground contamination, but in 
contrast to the IGM signal, the foregrounds are expe cted 



to be smooth in frequ ency, facilitating their removal ( Zal- 
darriaga et al.||2004 ). The 21 cm data will be supple- 
mented by further quasar absorptio n spectra (including 



clues from metal absorption lines: Oh 2002 Becker et al. 
20051), high redshift gamma ray bursts (e.g. Barkana 
& Loeb||2004 al |Totani et a l. 2005), high redshift galaxy 



surveys (e.g. [Kncib et al. 2004), and small-scale CMB 



measurements (Santos et al. 2003 Zahn ct al. 2005 Mc- 



Quinn et al.||20Tf5J), providing a wealth ol observational 



data on the process of reionization. 
Detailed theoretical modeling (see e.g . 



Gnedin 2000 



Razoumov et al . 2002; Ciardi et al. 2003; Sokasian et al.| 
20021 ISokasian et al.||2001MFurlanetto et al.||2004b|a| 



Kohler et al.|2005||iliev et al.|2005| IMellema et al.|2006 



is, however, required to constrain the topology of reion- 
ization from current and future observations. In partic- 
ular, it is subtle to infer quantities like the volume filling 
factor and size distribution of ionized regions, the correla- 
tion between ionized regions and large-scale over-density 
(i.e., does reionization proceed outside-in or inside-out), 
and the nature of the ionizing sources, from observations. 

Unfortunately, numerical modeling of reionization is 
challenging, requiring treatment of radiative transfer, 
preferably some treatment of gas dynamics, and a large 
dynamic range. A large dynamic range is required to 
resolve the small mass galaxies which may make up the 
sources and sinks of ionizing photons, while simultane- 
ously sampling the distribution of HII regions, which may 
be as large as R ~ 20 comoving Mpc hr 1 t owards the 
end of reionization ( |Furlanetto &i Oh 2005). For this 
reason, most calculations have been performed in pro- 
hibitively small simulation boxes, or been entirely ana- 
lytic (e.g. |Furlanetto et al.|2004a[ hereafter FZH04), al- 
though there has been very recen t progress towards large 
volume reio nization simulations ( |Kohler et al.|200"5} |Iliev| 
et al.|2005 1 . Indeed a skeptic might posit that, given the 
difficulty of simulating reionization, we will observe the 
21 cm signal before we can predict it. 

In this paper, we push forward by running a large vol- 
ume radiative transfer simulation. Our work represents 
progress on several fronts. First, we simulate reionization 
in a larger volu me tha n most previous works (although 



see Kohler et al. 2005; Ilicv et al.|2005 l, while maintain- 
ing niglTTn^ssTesoTuTfon - ^ThiTaTlows us to reliably calcu- 
late the size distribution of HII regions as well as power 
spectra of ionization and 21 cm fields, impossible with 
previous small volume simulations. Second, we compare 
our results with analytic calculations based on FZH04. 
These models are now widely used, and while elegant and 
inspir ed by previous small volum e reionization simula- 
tions ( Sokasian et al.|[2003 2004), they remain untested. 
Our comparison also gauges the level of theoretical con- 
trol in our modeling of reionization - i.e., how robust are 
our conclusions to the details of our modeling? One con- 
vincing way to dissuade the above-mentioned skeptic is 
to demonstrate that we can understand the gross features 

5 http://www.lofar.org 

6 http://www.skatelescope.org/ 



of our radiative transfer simulations analytically. Addi- 
tionally, if analytic models are sufficiently accurate then 
they are useful tools to forecast constraints from future 
experiments, and to construct mock surveys, providing 
important tests of data analysis procedures. This is im- 
portant given our ignorance of the nature of the ionizing 
sources: we would like to cover a large parameter space 
in the source properties, prohibitive with time-consuming 
radiative transfer simulations. Furthermore, future sur- 
veys will span volumes of several cubic Giga-parsecs, a 
challenging task for detailed simulations. 

We emphasize that our present work is only a first step 
towards more realistic simulations of Hydrogen reioniza- 
tion. As we describe subsequently, our radiative trans- 
fer simulations miss potentially important aspects of the 
physics of reionization. Specifically, we include only a 
crude prescription for the sources of ionizing photons, 
our coarse resolution underestimates the importance of 
recombinations - especially if mini-halos are present dur- 
ing reionization (|Haiman et al. 2000 Barkana & Loeb 
2002 Shapiro et al.|2Ulil ) - and misses small galaxies that 
may contribute ionizing photons, and we ignore feedback 
effects entirely. We intend to model some of these effects 
in the near future ( McQuinn et al.|[2006 )). 

Our wor k has overlap wit h the recent simulation and 
analysis of Iliev et al. (2005). In comparison to these au- 
thors, reionization finishes significantly later in our sim- 
ulation, near z ~ 6.5, as compared to z ~ 12, a con- 
sequence of our more conservative prescription for the 
ionizing sources. Moreover, our main present emphasis 
is in comparing our radiative transfer simulation results 
with 'hybrid simulations' based on analytic models. 

The layout of this paper is as follows. In SJ2]we describe 
our N-body simulation, source prescription and radia- 
tive transfer calculation. In ^3] we describe our 'analytic 
model simulation', which is more precisely an implemen- 
tation of a model based on FZH04 into the cosmological 
realization used for the radiative transfer simulation. We 
will sometimes refer to this scheme losely as an 'analytic 
calculation' although the implementation of the model is 
entirely numerical. In £J4]we present a detailed statistical 
description of our radiative transfer and analytic results. 
We describe a numerical scheme that incorporates the 
stochasticity of the source distribution into our analytic 
calculations in <|5J We also show that if extremely bright 
and rare sources reionize the IGM, bubble growth is less 
collective than in our fiducial model. 

In 36]we compare radiative transfer and analytic model 
predictions for the 21 cm signal. We conclude in SjTJ 
mentioning future research directions and emphasizing 
possible improvements to our simulations. 

Throughout we assume a fiat, ACDM cosmology pa- 
rameterized by: D, m — 0.3, Q,\ = 0.7, 0(, = 0.04, 
H = lOO/i km/s/Mpc with h — 0.7, and a scale-invariant 
primordial power spectrum with n = 1, normalized to 
( j 8 (z = 0) = 0.9[] 

7 This value for erg is slightly different tha n the value pr eferred 
by the WMAP satellite alone of 0.76 ± 0.05 ( |Spergel et al.|2006| ). 
When combined with other dat asets. such as the Lyman-a lorest 
| |Lewis|[2006| |Seljak et al.|[2006|, and weak lensing (see e.g. sec- 
tion 4.I.Y and table b ot jspergel et al.|2006| l , a higher value of <Tg 
can be found. Furthermore, changes in the fluctuation amplitude 
within the present experimental boundaries can be incorporated 
into our analysis by adjusting slightly the ionization efficiency pa- 



2. SIMULATIONS 

We begin by running a large N-body simulation to lo- 
cate dark matter halos, and produce a cosmological den- 
sity field. Next, we populate the dark matter halos with 
ionizing sources, usin g a simple prescription to connect 
mass and light (£2.2). In a subsequent post-processing 
step, we perform a radiative transfer calculation, cast- 
ing rays of ionizing photons fr om our sources through 
the cosmological density field ([2.3). We make two ap- 
proximations with this approach. First, we assume that 
the gas distribution perfectly traces the dark matter dis- 
tribution, as characterized by our N-body simulation. 
Second, we neglect the interplay between gas dynam- 
ics and radiation transport - i.e, in reality, structure 
formation responds to the passage of ionization fronts, 
and gas motions in turn influence the propagation of the 
fronts. These effects are essential in calculating the de- 
tailed small-scale behavior of ionization fr onts, as fronts 



slow down upon impacting dense clumps ( [Shapiro et al 
20041, but are less important for our goal of capturing 
the large-scale size distribution of HII regions. 

2.1. N-body simulations 

As noted in the introduction, we require a cosmo- 
logical simulation with a large dynamic range, in or- 
der to adequately sample the distribution of HII re- 
gions, while simultaneously resolving small galaxies. Ide- 
ally, we would resolve halos with virial temperatures of 
Tvir ^ 10 4 K - corresponding to a dark matter halo mass 
of M(j m ~ 10 8 M Q at z ~ 6 - above which atomic line 
cooling is efficient. In halos more massive than this, gas 
can cool, condense to form stars, and produce ionizing 
photons. This 'cooling mass' therefore represents a plau- 
sible guess as to the minimum host halo mass for ionizing 
sources. If molecular Hydrogen cooling is efficient de- 
spite radiative feedback, however, even smal ler mass ha- 
los should host sources (H aiman et al.fl997 |. Presently, 
we ignore this possibility. Additionally, high resolution is 
required to capture the dumpiness of the IGM, and prop- 
erly account for recombinations during reionization. On 
the other hand, HII regions may be l arger than R > 20 
Mpc Yc 1 at the end of reionization (Furlanetto & Oh 
2005| , necessitating a large volume simulation. Unlortu- 
nately, to resolve a 10 8 M Q halo with 32 particles, in a 
simulation box of side-length L = 100 Mpc h~ l , for ex- 
ample, requires a prohibitively large number of particles, 
N p ~ 3360 3 ! 

Our present N-body simulation is meant to represent 
a compromise between these competing requirements of 
large volume, and high mass resolution. Specifically, our 
N-body simulation follows 1024 3 dark matter particles 
in a box of side-length, L = 65.6 Mpc h~ l , usi ng an en- 



hance d version of the TreePM code, Gadget-2 ( [Springel 
2005 1 . We run the simulation assuming the flat LCD1V1 



cosmology specified in the introduction, with initi al con- 
ditions generated using the Eisenstein & Hu ( 1999 1 trans- 
fer function. 

Dark matter halos are identified from simulation snap- 
shots, usin g a friends-of- friends algorithm (e.g., Davis 
et al.|1985 1 . Specifically, particles are grouped into halos 
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Fig. 1. — Halo mass function from our N-body simulation. The 
black points with (Poisson) error bars indicate the halo mass func- 
tion from our simulation as a function of redshift. The green curve 
is the Sheth-Tormen fitting function for the halo mass function, 
while the red dashed line shows the Press-Schechter fitting func- 
tion. 



using a linking length of b = 0.2 times the mean inter- 
particle separation. Linked groups of greater than 32 
particles are considered resolved, and to constitute dark 
matter halos. This corresponds to a minimum halo mass 
of 1O 9 M , just an order of magnitude above the cooling 
mass. 

The resulting mass function is shown in Figur e [Tj span- 
ning a b r oad red shift range between z ~ 6 — 20. |Barkana] 
& Loeb ( 2004b ) showed that if a simulation is normal- 
ized to the cosmic mean density, the ha lo mass functions 
will be biased. According to Figure 3 oflBarkana fc Loeb| 
(2004b), the bias introduced in our calculations should 
only be of order 0.1 %. The halo mass function is sam- 
pled with large dynamic range, roughly three orders of 
magnitude near z ~ 6. The simulated mass function 
is always l arger than predicted by t he Press-Schechter 



formalism (Press & Schechter 1974), bu t generally in 
good agre ement with the Sheth-Tormen ( ]Sheth fc Tor 
men|1999] ) fitting formula. At the highest redshitts sam 
pled, however, our results fall in between the two fitting 
formula. This is in qualitative agreem ent with recent 



measuremen ts from Reed et al. (2005[), and Heitmann 



rameter. This does not qualitatively affect our results, as we con- 
firmed within the analytic scenario. 



et al. (2006), although our ma ss function appea rs sys- 
tematically higher than that of Iliev et al. ( 2005 ) . The 
figure shows that the abundance of our lowest mass halos 
is systematically below theoretical expectations, likely a 
consequence of our limited mass resolution. As a conser- 
vative measure, we therefore place ionizing sources only 
in halos of mass larger than M m ; n = 2 x 1O 9 M0, corre- 
sponding to a 64— particle halo. 

2.2. Ionizing Sources 



Our next step is to connect mass with light that is, 
we wish to populate the dark matter halos from our N- 
body simulation with ionizing sources. In this paper, 
we will adopt a very crude prescription for our ioniz- 
ing sources, leaving a more sophisticated prescription to 
future work. This will facilitate comparison with the an- 
alytic models (see jJ3J). Specifically, we populate each 
dark matter halo with a single source whose luminos- 
ity in Hydrogen ionizing photons is directly proportional 
to the host halo mass, N = cMhaXo- Clearly the pa- 
rameter c encodes a good deal of complicated physics, 
involving the efficiency of star formation, the efficiency 
of producing ionizing photons, the fraction of ionizing 
photons that escape from the host halo, etc. With this 
single simplifying assumption, the cumulative number of 
ionizing photons released by the sources, per hydrogen 
atom in the IGM, at time t is Np^/Nn oc f Q dt 1 f co \\(t'). 
Here /coii(t') is the fraction of mass i n halos with mass 
M > M min = 2x 10 9 M o . Using the |Sheth fc Tormen 



( 1999 ) mass function, which closely matches our simu 
lation results (Figure [l]), we find that c = 3.1 x 10 41 
photons/sec/M Q yields onephoton per hydrogen atom 
at z = 6.5. (See Figure [2] and associated text for a 
discussion). This choice of c corresponds roughly, for 
example, to Pop II stars, forming with an efficiency of 
/* = 0.1 from a Salpeter IMF, with a stellar lifetime 
of At ~ 5 x 10 7 yrs, and a modest escape fraction of 



/esc ~ 0.01 ( Loeb et aL]|2005 1 . We adopt this conversion 
in all subsequent calculations. 

2.3. Radiative Transfer 

We next form a coarse density field for many snapshots, 
spaced in equal time intervals of At — 5 x 10 7 years and 
spanning a broad redshift range from z ~ 6— 16, by grid- 
ding our dark matter particles onto a unifo rm. Cart esian 
grid with 256 3 mesh points. Our sources (S2.1 and S2.2| 
are tabulated at the same time-sampling. 



and 



moved 

close to the center of their corresponding cell. Occa- 
sionally, several sources land in a single cell and our con- 
sidered to be a sole, more luminous source. At z ~ 6.5, 
near our assumed completion of reionization, there are 
~ 330, 000 ionizing sources in our simulation. 

With the ionizing sources and cosmological density 
field in hand, we trace rays of ionizing photons through 
the simul ation box using th e adaptive ray-tracing scheme 
of |Abel h Wa ndelt (20Q2 J), and the code of So kasian et 
al. ( |Sokasian et al.pOTJTnSokasian et al |2003[). Further 
impro vements to the code were made (McOuinn et al. 
20061. We refer the reader to these papers for the de- 
tails of the code used here, but give a brief summary. 
In short, the code assumes a sharp ionization front, and 
tracks the position of the front by casting rays a nd in- 
tegrating ov er the ionization front jump condition ( Abel 



et al. 19981. The jump condition amounts to tabulat 
ing the number of photoionizations and recombinations 
along a ray, halting the ray when its photon supply is ex- 
hausted. Each source is considered separately, although 
the order in which sources are pr ocessed is randomize d 
at each timestep to avoid artifacts ( Sokasian et al.|2003 ). 
Behind the ionization front, each source hitting a 
given cell contributes a photoionization rate of Thls = 
aN /(4irr^), i.e. assuming optically thin conditions 
within the front. Here r s is the distance from the cell 



3 0.5 




Fig. 2. — Ionization fraction as a function of redshift. The black 
circles show the mass-weighted ionization fraction from the simu- 
lation, while the blue squares show the volume-weighted ionization 
fraction. The red line is the cumulative number of ionizing pho- 
tons per hydrogen atom expected for our ionizing sources. The 
close resemblance between the number of photons per atom and 
the measured ionization fractions owes to the poor resolution of 
our radiative transfer calculation, which underestimates the im- 
portance of recombinations. 



in question to a source, N is the number of Hydrogen 
photons per second from a source, and a is a frequency- 
averaged cross section, comput ed here assuming each 
source has a spectrum ex v~~ A ( |Sokasian et al. 2001). 
Within the front, ionization fractions are computed as- 
suming ionization equilibrium and a uniform tempera- 
ture of T = lO^K, and neglecting sub-grid c lumping. 
We follow the approach of ( Sokasian et al.|20"oT ) in using 
case B recombination rates when casting rays through 
the grid, while using case A coefficients to compute the 
ionization fraction in previously ionized cells. Helium is 
assumed to be at most singly-ionized by our soft sources, 
and we assume that the Hell front precisely tracks the 
HII front. Similarly, inside the front we assume that the 
ionized Helium (Hel l) fraction traces the HII fraction 
( |Sokasian et al.|2003 ). Note that all of these assumptions 
impact mainly the detailed ionization fractions within 
the front, and are less important for tracking the overall 
size d i stribu t ion of HII regions. In contrast to Sokasian] 
et al. (2001); Sokasian et al. (2003), we do not include a 
diffuse background radiation field, simply allowing rays 
to wrap around the periodic box. 

Our assumption of a sharp ionizing front is justified 
given the short mean free path of Hydrogen ionizing pho- 
tons in the pre-reionization IGM. Mellema et al. (2005) 
present explicit comparisons between 'ionization front 
tracking' and more detailed calculations that self con- 
sistently solve for the optical depth, ionization fraction, 
and temperature. At least in the case of a single source 
(their Figure 16), ionization front tracking reproduces 
very closely the results of more detailed calculations, fur- 
ther justifying our approach. 

In Figure [2] we plot the redshift evolution of the ioniza- 
tion fraction in our simulation. The black circles show 



the mass-weighted ionization fraction, while the blue 
squares show the volume-weighted ionization fraction. 
The mass- weighted ionization fraction is somewhat larger 
than the volume- weighted ionization fraction. This is be- 
cause the ionizing sources in our simulation are highly bi- 
ased, and ionize their overdense envi rons before breaking- 



free to ionize nei ghboring voids (e.g. Sokasian et al. 2003 t 
Iliev et aT]|2005 1 . The reionization process takes a fairly 



significant stretch of cosmic time, with the mass- weighted 
ionization fraction at the level of x-^ m ~ 0.1 at z ~ 9, and 
attaining x lm ~ 1 only by z ~ 6.5. 

The evolution of the neutral fraction in our model is 
consistent wi th, although not required by the measure- 
ments of e.g. Fan et al. (2005), which demand only that 
the IGM reiomze sometime before z > 6. Our model pro- 
duces an electron scattering optical depth of r e = 0.06 , 
on the low side of CMB constraints (Page et al. 2006), 
which sugge st t„ = 0.09 ± 0.03. We emphasize that our 



choice of c ({ 2.2 1 was calibrated so that reionization ends 
slightly above z > 6, so this should be viewed as a con- 
sequence of our assumptions, rather than a theoretical 
prediction. Although our model is tuned to give late 
reionization, analytic models find that the size distribu- 
tion of HII regions depends primarily on the bias of the 
ionizing sources, with only an implicit dependence on 
redshift (Furlanetto et al. 2006a). The size distribution of 



HII regions at a given ionization fraction is therefore ex- 
pected to be a robust result, independent of our detailed 
assumptions about the efficiency of th e ionizing sources, 
(although see Furlanetto et al. 2006a, fjslfor caveats). 

Note that our simulation terminates slightly before 
reionization completes (xi(z) ~ 1). We stop our calcu- 
lation early because we do not include a 'diffuse back- 
ground' in our simulation (Sokasian et al. 20021, and 
so our calculation becomes very expensive at the end 
of reionization when rays wrap around the simulation 
box several times. In any event, the tail end of reion- 
ization is likely poorly modeled in our simulation, since 
this stag e may be regulated primarily by Lyman limit 
syste ms ( Miralda-Escude et al.||2000" Furlanetto & Oh 
20051, which are missing in our analysis. 

In our simulation, the mass-weighted ionization frac- 
tion closely tracks the cumulative number of ionizing 
photons per Hydrogen atom emitted by our ionizing 
sources (see the solid red line in Figure J2J), but this is 
partly an artifact of the poor resolution of our radiative 
transfer calculation. Our low grid resolution underesti- 
mates the amount of small scale structure in the density 
field, and hence the importance of recombinations. In 
the future, we intend to model recombinations as 'sub- 
grid physics', accounting for enhancements in the recom- 
bination rate_ow2n^tojiiiresolved small scale structure 
(see e.g. |Kohler et al. (2005 )). Presently, we caution that 
we are under-estimating the number of ionizing photons 
per Hydrogen atom required to complete reionization. 
Furthermore, we expect reionization to be even more 
extended than in our calculation, since recombinations 
should slow the growth of HII regions. 

3. NUMERICAL SCHEME BASED ON ANALYTIC 
CONSIDERATIONS 

As motivated in the introduction, we compare our re- 
sults with a hybrid scheme inspired by the analytic model 
of FZH04. In this section, we describe this hybrid model, 



but also refer the reader to FZH04, Zahn et al. (2005 1 for 
more information. The major advantage of our imple- 
mentation over a purely analytic calculation is that the 
hybrid scheme, which amounts to a Monte-Carlo realiza- 
tion of the analytic model, can capture the asphericity 
of HII regions during reionization. 

The hybrid scheme starts by considering spheres of 
varying radius surrounding every point in the IGM. 
Within each such sphere, we calculate the total ioniz- 
ing photon yield of the sources, and the total enclosed 
mass in neutral Hydrogen. In the event that the photon 
yield in a sphere around a given point exceeds the num- 
ber of interior Hydrogen atoms, the point is considered 
ionized. FZH04 show that this amounts to a barrier- 
crossing proble m, solvable with to ols from the excursion 



set formalism (Bond et al. 1991). Specifically, we as 



sume that the mass contained in halos in a region of total 
mass m, and over-density <5 m , follows the extend ed Press 



Schechter formula for the collapse fraction (e.g. Lacey & 
CoIelp93l ): 



/coll. (m > m min \S m , z) = erfc 



5 c (z) - S n 



V2[al in -aHm)} 



Here a 2 (m) is the (present day) linear variance of density 
fluctuations on the scale m, 5 c (z) = 1.686/ D(z) is the 
critical density for collapse scaled to today, and D(z) is 
the linear growth factor. The quantity <jf nin is the linear 
variance smoothed on a mass scale corresponding to that 
of the minimum mass halo that can host ionizing sources, 
presently m mln = 2x 10 9 M Q (jJ2T|. 

For constant mass to light sources, the criterion for a 
region to self-ionize is then: 



dt' /coil. (m > m m in\S m ,t') > 1 , 



(2) 



where a is an efficiency factor linking halo mass and ion- 
izing photon yield. A region can self-ionize if it is suf- 
ficiently overdense to satisfy the inequality in Equation 
(pi) . Note that this is a slight modification from FZH04 to 
the case of ionizing sources with a constant mass to light 
ratio, as assumed in our simulation. In practice, however, 
we find that the threshold criterion of Equation (pi) gives 
quantitatively similar results to that of FZH04, although 
it produces slightly larger HII regions. 

Our hybrid scheme then amounts to smoothing the lin- 
ear density field generated from the initial conditions of 
our N-body simulation, and checking whether cells sat- 
isfy the condition of Equation (pi). Algorithmically, we 
start by considering large spheres( comparable to the size 
of our simulation box) , and gradually stepping down in 
radius, eventually reaching smoothing scales compara- 
ble to that of our simulation pixels. At each radius we 
keep track of which cells satisfy the condition of Equa- 
tion (pi). In the event that a cell does not cross the bar- 
rier, i.e. satisfy the condition of Equation (pi), at any 
smoothing scale, the cell is considered neutral. Proceed- 
ing from large smoothing scales and progressing down- 
ward to smaller smoothing scales ensures that we account 
correctly for cells that are ionized by neighboring sources 
(FZH04). At this stage, we have a map of the ionization 
field consisting of T's (completely ionized pixels), and '0's 



G 



(completely neutral pixels). This is a good approxima- 
tion to the true equilibrium ionization fractions, which 
are expected to be very close to unity. 

This scheme is quite fast: for our present 256 3 grid cal- 
culation, with 50 logarithmic smoothing steps, the com- 
putation (at a given redshift) takes only ~ 12 minutes 
on a desktop computer with a 3 GHz processor. This is 
vastly more efficient than our full radiative transfer cal- 
culation: our N-body simulation takes 38 hours to run 
down to z ~ 6 using 134 2 GHz processors, and our post- 
processing calculation requires a few additional days of 
running time on a large memory computer. With our 
rapid numerical scheme, we can produce an ionization 
map based on the analytic model and compare with our 
radiative transfer simulations. Using precisely the ini- 
tial conditions from our N-body simulation in our hybrid 
calculation allows us to compare radiative transfer and 
analytic ionization fields on a cell-by-cell basis. 

Before presenting this comparison, there are a few 
more pertinent technical details. Ideally, we would com- 
pare the analytic and radiative transfer calculations with 
identical assumptions regarding the ionizing efficiency of 
our sources, i.e. we should calibrate a in Equation (f2j) 
based on the source prescription of |2.2| In practice there 
are several difficulties with matching precisely the sim- 
ulated source prescription. Most important, Equation 
(nj) is derived assuming sharp fc-space filtering, while our 
smoothing procedure adopts a spherical top-hat in real 
space. This slight inconsistency in our modeling means 
that our model does not conserve photons precisely, af- 
fecting the ionization fraction for a given source effi- 
ciency, a (see the Appendix). Further, our simulated 
mass function is closer to the Sheth-Tormen fitting for- 
mula (Shcth & Tormcn 1999) than the Press-Schechter 
(Press & Schcchter] |l974| mass function, and we re- 
quire an analo gue of Equation ([!]) for the Sheth-Tormen 
mass function ( Barkana fc Loeb|f2004b Furlanetto et al 

P). We improve on some of these shortcomings in 
le upshot of this is that, in order to compare with 
our radiative transfer simulations, we adjust a in Equa- 
tion (p) at each redshift to match the (volume-weighted) 
ionization fraction. This readjustment is usually of order 
20%. 

We show examples of the resulting ionization maps 
in Figure [3J The left column shows thin slices through 
the radiative transfer simulation at three different stages 
in the reionization process: z = 8.16, 7.26 and 6.89 
when the (volume-weighted) ionization fraction is x- ltV = 
0.11,0.33, and 0.52 respectively. The right column shows 
corresponding slices from the hybrid simulation scheme. 
Several conclusions are immediately apparent. 

First, the ionized regions are quite large at the inter- 
mediate and late stages of reionization. The ionizing 
sources are highly clustered, and HII regions quickly start 
growing collectively around the sources, rapidly reach- 
ing much larger sizes than can be achieved by individ- 
ual sources (FZH04) r] Second, the hybrid simulation 
is in good general agreement with the radiative trans- 
fer simulation. The hybrid scheme seems to 'locate' the 



8 The size of an HII region belonging to an individual average 
source would be too small to display in this Figure, e.g. for a 10 
Mq source shining for one simulation time step of 5 • 10 7 years this 
size would be only 0.3 Mpc/h. 
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Fig. 3. — Maps of the ionization field. The left column shows HII 
regions for a thin slice through our radiative transfer simulations 
at redshifts z=8.16, z=7.26, and z=6.89 (top to bottom). The 
volume-weighted ionization fraction at these redshifts is Xi v = 
0.11, 0.33 and 0.52, respectively. The slices are 0.25 Mpc//i deep, 
and 65.6 Mpc/h on a side. The right panel shows the same using 
our hybrid simulation scheme, as applied to the initial conditions 
used in our radiative transfer simulation. The analytic modeling 
agrees well with the more detailed simulation, although there is 
more small scale structure in the map from the radiative transfer 
simulation (see text). 



HII regions found in the radiative transfer calculation, 
and additionally reproduces their general morphology. 
Third, the HII regions in the analytic calculation are a 
bit more 'connected' than those in the radiative trans- 
fer simulation. Equivalently, the ionization field in the 
radiative transfer simulation appears to have more small 
scale structure than the ionization field from the hybrid 
scheme. In the following sections, we will quantify the vi- 
sual comparison of Figure [3] diagnose differences found, 
and refine our numerical scheme. We contrast the mor- 



phology seen here with that from Iliev et al. (2005) in 
4. STATISTICAL DESCRIPTION 



In this section we present a detailed statistical descrip- 
tion of our results. Throughout we will compare with 
our hybrid scheme rather than the purely analytic calcu- 
lations for two reasons. First, there are technical difficul- 
ties in the analytic calculations at intermediate ionization 
fractions (McQuinn ct al. 2005), and second, we would 
like to be able to model non-spherical bubble shapes. 

4.1. The Bubble PDF 

The first statistic we consider is the probability distri- 
bution of bubble sizes. That is, we calculate how large 
the HII regions are at different stages of reionization. 
This depends somewhat on how one chooses to define 
contiguous ionized volumes - Figure [3] clearly illustrates 
that the ionized regions are not spherical, particularly 
at the end of reionization. The ionized regions do, how- 
ever, obtain a reasonably well-defined characteristic size 
at each redshift. In order to quantify this, we require a 
convenient and well-motivated definition of 'bubble' that 
we can apply consistently to the radiative transfer simu- 
lation and the hybrid scheme. 

Here we adopt a definition of bubble size inspired by 
the excursion set formalism, upon which our analytic cal- 



culation is based (see Iliev et al 



2005 for an alternate 



approach). Specifically, we 'draw spheres around each 
point in our simulation box of varying radius, i?, and 
average (smooth) the ionization field within each such 
sphere. We start by considering large spheres, of vol- 
ume comparable to that of our simulation box, and step 
downward in size until we eventually get to the size of our 
simulation pixels. At each smoothing radius, R, we com- 
pare the average ionization in each sphere to a threshold 
ionization, x t h- A pixel is marked as 'ionized' and be- 
longing to a bubble of radius R, when R is the largest 
smoothing radius at which the pixel's smoothed ioniza- 
tion exceeds the threshold ionization, x t h- If a given pixel 
fails to exceed the threshold ionization at all smoothing 
scales, it is considered neutral (not ionized). 

The bubble pdf is then derived by tabulating the frac- 
tion of ionized pixels that lie within bubbles with radius 
between R and R+dR. With this convention, the bubble 
pdf is normalized to unity rather than to the mean ion- 
ization fraction. The results of this calculation are shown 
in Figure|4j for an ionization threshold of 2th = 0.9. The 
figure illustrates quantitatively the visual impression of 
Figure [3] : the HII regions have a well-defined character- 
istic size at each stage of reionization, and this character- 
istic scale evolves as bubbles around ne ighboring sources 
overlap and grow collectively (FZH04, Furlanetto et al. 
2006a). The characteristic scale evolves from sub-Mpc 
scales at z = 8.16, when the volume- weighted ioniza- 
tion fraction is xi iV = 0.11 to R > 10 Mpc co-moving 
at z — 6.56 when the volume- weighted ionization frac- 
tion is X[ yV = 0.77. The large size of HII regions at high 
ionization fraction implies that large volume simulations 
are r equired to adequately sa mple this stage of reioniza - 
tion ( |Barkana fc Loeb||2004b , FZH04, [iliev ' et al.l)2005| ). 
The precise value of the characteristic bubble size de- 
pends somewhat on the number we adopt for the thresh- 
old ionization. For instance, if we instead adopt the 
less stringent threshold of ccth = 0.7, the characteristic 
size increases by a factor of ~ 2 near z = 8.16. Again, 
while our definition of bubble-size is somewhat arbitrary, 
the bubbles nevertheless have a well-defined characteris- 
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Fig. 4. — Size distribution of HII regions as a function of red- 
shift. The solid curves show results from the radiative transfer 
simulation, while the dotted curves are from the analytic calcu- 
lation. We adopt a threshold ionization of x t ^ = 0.9 (see text). 
The volume-weighted ionization fractions at the redshifts shown 
are x ijV =0.11, 0.20, 0.33, 0.52, 0.77 at z = 8.16, 7.68, 7.26, 6.89 and 
z = 6.56 respectively. 



tic scale ( Furlanetto et al.|200 6a), and our algorithm can 
be applied consistently to each of the analytic model and 
radiative transfer ionization maps. 

The dotted lines indicate that our hybrid scheme repro- 
duces the bubble pdf simulated through radiative trans- 
fer quite accurately, roughly matching the characteris- 
tic bubble size and its trend with redshift. The hybrid 
scheme however leads to slightly larger HII regions at all 
but the final redshift. We will discuss this difference in 
future sections. At the final redshift, the agreement is 
almost exact, however here our simulated volume is too 
small to provide a representative sample. 

4.2. Power Spectra of the ionized fraction 

For further comparison, we measure the (spherically 
averaged) 3d ionization power spectrum as a function of 
redshift. We consider the ionization field S x — x(f) — (x), 
where x(r) denotes the ionization at spatial position r, 
and (x) denotes the volume-averaged ionization. Note 
that we do not normalize by the mean ionization here, i.e. 
we consider the absolute ionization fluctuation, rather 
than the fractional fluctuation. The result of the power 
spectrum calculation is shown in Figure [5] with power 
spectra calculated from the radiative transfer simulation 
plotted in red. Throughout this paper we plot the dimcn- 
sionless power spectrum, A 2 (/c) = k 3 P(k)/(2ir 2 ), which 
yields the contribution to the variance per logarithmic 
interval in fc. On large scales at high redshift the ion- 
ization power spectrum is proportional to the density 
power spectrum, while it turns over or flattens on scales 
in which there are ionized bubbles. On intermediate 
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Fig. 5. — Power spectra of the ionized fraction, going from large 
redshift (small ionization fraction) to small redshift (large ioniza- 
tion fraction). The red lines are from the radiative transfer simu- 
lation, the blue dashed lines are from the analytic hybrid calcula- 
tion, while the purple dotted lines show results from the improved 
scheme of the next section. The high-fc behavior (k > lOh Mpc -1 ) 
is an artifact from discreteness noise. 



scales the power spectrum from the radiative transfer 
simulation has a somewhat larger amplitude. We at- 
tribute this to a superior tracking of the density field 
around the edges of the bubbles: rays can travel into 
underdense regions and will be hindered by overdensi- 
ties. In the numeric schemes the features of the density 
field are washed out somewhat, resulting into a generally 
smoother edge structure, see also Figure [3J The bub- 
ble 'feature' moves to progressively larger scales (small 
k) as reionization proceeds, a further illustration of the 
bubble growth seen in Figure [I] The blue dashed curves 
show power spectra from ournybrid simulation, which 
are similar to the radiative transfer power spectra, ex- 
cept with slightly more large scale power, and slightly 
less small scale power. One can also infer from the figure 
that an even larger volume simulation is preferable, in 
order to better sample the large scale ionization power 
spectrum. Finally, the purple dotted lines are from an 
improved numerical scheme which we discuss in the next 
section. 

In order to further quantify the agreement between the 
radiative transfer simulation and the hybrid scheme, we 
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Fig. 6. — Cross correlation coefficient between the ionization 
fields from the radiative transfer simulation and the analytic model 
calculations. The thin lines show the cross correlation coefficient 
between the radiative transfer and hybrid simulations at a few dif- 
ferent redshifts. The thick lines show corresponding results from 
the improved hybrid simulation described in the next section. 



calculate the cross correlation coefficient between the two 
ionization fields. The cross correlation coefficient is de- 
fined by r{k) = A* liX2 (fc)/ [A^A^fc)] 1/2 . In this 
equation A xl x2 (fc) is the cross power spectrum between 
the radiative transfer simulation and hybrid scheme ion- 
ization fields, while A xl (fc), and A^ 2 (k) are their respec- 
tive power spectra. The cross correlation coefficient is 
bounded between 1 and —1, with r(k) = 1 indicating 
perfectly correlated modes, and r(k) — — 1 designating 
perfectly anti-correlated modes. The results of this cal- 
culation are shown as thin lines in Figure p] (ignore, for 
now, the thick lines which show results from the im- 
proved hybrid scheme introduced in the next section). 
The correlation coefficient is always larger than r ^ 0.5 
for scales larger than k < lh Mpc -1 , while it drops off 
on smaller scales. This quantifies the qualitative agree- 
ment suggested by Figure |3j the radiative transfer and 
hybrid scheme ionization fields trace each other closely 
on scales larger than k < lh Mpc -1 . The cross corre- 
lation between the two fields becomes slightly weaker at 
low redshift, as the average ionization increase. A plau- 
sible explanation for the slightly worse agreement at low 
redshift is that our hybrid simulation scheme has diffi- 
culty with 'bubble mergers' (see the Appendix), which 
are more frequent at high ionization fraction. 

Why does the cross correlation between the two fields 
drop off around k > lh Mpc -1 ? The analytic model 
assumes a one-to-one correspondence between the abun- 
dance of halos and the (Lagrangian) matter overdensity 
on a given smoothing scale. We know this is inexact. 
For one, the abundance of our minimum mass sources is 
Mdn/dM < 1 Mpc -3 . On ~ 1 Mpc scales, we therefore 
expect significant Poisson scatter in the abundance of 
ioni zing sources in our radiative transfer simulat ion (see 
also |Furlanetto et al.|2006a||Cohn fc Chang|2006l ) . To ex- 
plore this further, we compute the cross power spectrum 
between the halo density field and the matter density 
field. The cross correlation coefficient between the halo 
and matter density fields qualitatively mirrors the cross 



correlation between the two ionization fields seen in Fig- 
ure pi dropping off at fc > \h Mpc -1 . In other words, 
the halo bias is stochastic on scales of k > \h Mpc -1 for 
our assumed source population. This stochasticity is not 
incorporated in our analytic hybrid scheme, and likely 
leads to the lack of small scale structure compared to 
the ionization held simulated through radiative transfer. 
We will return to this issue in Sj5] We note here, how- 
ever, that this Poisson scatter would presumably be less 
important if our radiative transfer simulation resolved 
smaller, more abundant galaxies. 

The analytic model connects ionized regions with large 
scale overdensities, which contain more sour ces and are 
reionized bef ore underdense regions (FZH04, Barkana & 
Loeb 2004b). The model therefore predicts that the 
ionization held is positively correlated with the matter 
density (e.g. McQuinn et al. 2005b), before turning 
over on scales comparable to that of the ionized bub- 
bles (FZH04). Figure [7] shows the cross power spec- 
trum between ionization and density (bottom panel) as 
well as the cross-correlation coefficient between the two 
fields. The radiative transfer simulation results (solid 
lines) nicely mirror the analytic model predictions (dot- 
ted lines). Since the analytic model ionization field is 
based on the initial condition density field, it is slightly 
less correlated with the evolved density field than the ra- 
diative transfer simulation ionization field. In our radia- 
tive transfer simulation and hybrid scheme, reionization 
proceeds inside-out with the overdense regions reionized 
before underdense regions, as emphasized by FZH04 and 
Sokasian et al. (2003, 2004). Recombinations, under- 
estimated in our present simulations, could potentially 
weaken this correlation or, in an extre me case, reverse 
the correla tion with voids ionized first ([Miralda-Escude 
et al.|[2000). We intend to explore this in future work. 



5. IMPROVED NUMERICAL SCHEME 

Although the agreement between our radiative trans- 
fer simulation and the hybrid scheme is already quite 
good, we present here a modified numerical scheme that 
impro ves upon the one presented in S|3] and |Zahn et al.| 
(2005). Specifically, we aim to fix two short-comings 
of the analytic calculation. First, as mentioned previ- 
ously, the analytic calculation is based on the Press- 
Schechter formula for the collapse fraction. This for- 
mula is derived assuming sharp fc-space filtering, while 
our scheme filters the initial density field with a top- hat 
in real spa ce, which is slightly inconsistent (McQuinn 
et al. 2005 the Appendix). Second, the mass function 
m our radiative transfer simula tion (Figure nj is closer 
to the Sheth & Tormen (1999) mass function than the 
Press & Schcchtcr (1974| ) mass function. Finally, the an- 
alytic calculation assumes a one-to-one correspondence 
between initial over-density and halo abundance. As we 
discussed in j ]4.2| the halo bias in our N-body simulation 
is stochastic on small scales. 

Each of these shortcomings can be remedied by directly 
using the simulated halos in our numerical scheme, rather 
than the Press-Schechter formula for the collapse frac- 
tion. More specifically, we place the halo distribution 
from our N-body simulation on a grid and compare, at 
each grid cell, the halo mass to the total mass enclosed 
by a spherical top-hat. We then use a condition analo- 
gous to Equation (pi) to determine whether a region is 
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Fig. 7. — Top panel: Cross correlation coefficient between the 
ionization and density field. The solid (dotted) lines show the cross 
correlation coefficient between the ionization and density fields in 
the radiative transfer simulation (hybrid scheme) at several red- 
shifts. Bottom panel: Cross-power spectrum between the ioniza- 
tion and density field. Solid lines are calculations from the radiative 
transfer simulation, while dotted lines are from the hybrid scheme. 



ionized by the sources within it. In other words, the cal- 
culation proceeds exactly as in f|3j except that we use 
the halo distribution directly from the simulation, rather 
than Press-Schechter theory. Note further that we now 
consider the evolved, non-linear density field rather than 
the initial, linear density field to determine if a region can 
self-ionize. We will call this improved numerical imple- 
mentation the 'halo-smoothing' scheme in what follows. 
The CPU intensity of this scheme is again dominated by 
the number of FFT's necessary to achieve convergence 
in the bubble size statistic. As with the analytic scheme, 
this is roughly 12 minutes on a 3GHz Intel Xeon desktop 
computer. 

The results of this new scheme are shown in compar- 
ison w ith radiative transfer and analytic calculation in 
Figure [TO} where we show 21 cm brightness temperature 
fluctuations (see ^6| for a thin slice through the simula- 
tion volume. This is analogous to Figure [31 except the 
ionized regions are now dark, the neutralregions now 
bright, and fluctuations in the gas density are now visi- 
ble in the neutral regions. The left column shows results 
from our radiative transfer simulation, the right column 
shows the standard FZH04-type implementation, while 
the center column shows our improved halo-smoothing 
scheme. The blue dots in the left and center column 
show the ionizing sources contained in the thin simula- 
tion slice. The new scheme clearly resembles the full 
simulation more closely, with more disconnected ionized 
regions, owing to the presence of Poisson fluctuations in 
the source distribution. 

Figure [5] quantitatively illustrates improved agreement 
with the radiative transfer calculation, with our im- 
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proved scheme showing more small scale power than the 
hybrid simulation scheme. Figure [6] additionally shows 
the cross correlation between the radiative transfer ion- 
ization field and the ionization field in the improved nu- 
merical scheme (thick lines) . The halo-smoothing ioniza- 
tion field traces the ionization field simulated through ra- 
diative transfer more closely, and down to smaller scales, 
than in our initial calculation. We attribute the improved 
agreement largely to our incorporation, in the improved 
scheme, of Poisson scatter in the halo abundance. 

If the ionizing sources are even less abundant than we 
assume presently, the Poisson scatter naturally becomes 
more important. Indeed for sufficiently rare sources, 
Poisson fluctuations dominate over source clustering on 
the scale of a typical bubble, and bubble growth is less 
'collective' than in our fiducial model. In this regime, 
the morphology of HII regions during reionization may 
be qualitatively different. To examine this, we repeat our 
halo-smoothing calculation at z = 7.26 including only ha- 
los with m > 4x 10 10 M Q as sources. We adjust the ioniz- 
ing efficiency of these rarer sources upward to match our 
usual ionized fraction at this redshift, x- uv = 0.33, in or- 
der to compare maps at fixed ionization fraction. The re- 
sult of this calculation is shown in the left panel of Figure 
One can see that the bubbles are considerably more 
spherical than in our usual source prescription (middle 
panel), and that the HII regions have a more sharply 
defined scale. The left panel further illustrates that for 
this source prescription there are very few sources in each 
bubble. Note that this is a thin slice, and some sources 
contributing to bubble growth lie above or below it. 

Furthermore, the left panel qualitatively resemble s the 
morphology seen in the reionization simulations of Iliev 
et al. ( 2005 ) (see their Figure £rl) ■ Their simulations are 
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done at higher redshift, but have a similar source number 
density as our present, extreme choice of m > 4x 10 10 Af Q 
(with this choice our simulation volume contains roughly 
5, 000 sources a t z=7.2 6). We regard the morphology seen 
in Iliev et al. ( 2005 1 as unlikely to represent the true 
morphology of Hll regions during reionization. Their 
choice of minimum source mass (M m [ n = 2.5 x 1O 9 M ) 
is driven by the low mass resolution of their simulations, 
and the efficiency of their ionizing sources is boosted ex- 
tremely high in order to m atch first-year WMAP con- 
straints (Kogut et al. 20031. In other words, their sim- 



ulation represents a very extreme case of reionization 
by rare, bright sources. Our simulation is also missing 
plausible ionizing sources, given our comparable mini- 
mum source mass. However, owing to the different as- 
sumptions about the ionizing efficiency in our simulation, 
reionization occurs later and so our sources are much 
more abundant (265,000 sources in the simulation vol- 
ume at z=7.26). We are hence still in the regime where 
HII regions grow collectively, and we expect only small 
modifications to the morphology and size distribution of 
HII regions when we include still smaller mass sources. 
This is illustrated in the right panel of Figure [8] where we 
show predictions for our original hybrid scheme (§31), with 
the minimum source mass extended down to the cooling 
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While there are some differ- 



9 Note for comparison with these authors' figure: their simu- 
lation box has a side length of L = 100 Mpc//i, while ours has 



ences with the results from our usual source prescription 
(center panel) , the differences are clearly smaller than in 
comparison to the Pois son-dominated cas e (left panel). 
^\ The differences with [Iliev et al.| ( |2005| highlight the 
utility of our fast numerical schemes tor quickly examin- 
ing many different prescriptions for the ionizing sources 
and for understanding the robustness of the results. 

6. 21 CM SIGNAL AND POWER SPECTRA 

The statistics discussed in £|4] are largely diagnostic, 
aimed at describing the size distribution of HII regions 
in the simulation, and characterizing the agreement be- 
tween the radiative transfer simulation and analytic cal- 
culations. In this section we make a more observationally 
relevant comparison, contrasting radiative transfer and 
analytic 21 cm power spectra. 

The 21 cm brightness temperature, relative to the 
CMB, at observed freq uency, v, and redshift, z, is (e.g. 
Zaldarriaga et al.||2004[ ): 
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where S s is the density contrast of gas in redshift space, 
and Ts is the spin temperature of neutral hydrogen. At 
the redshifts we consider presently, the 21 cm excita- 
tion temperature is likely coupled to the gas temperature, 
and much larger than the temperature of the CMB (e. g. 
Furlanetto |2006"1 |Chen fc Miralda-Escude|2004||Ciardi 



Madau||2003l 7Ts » r CM B, implying 6T oc (1 + 6 s )x H 
We then model the 21 cm brightness temperature us- 
ing the simulated density and peculiar velocity fields, in 
conjunction with radiative transfer/analytic calculation 
simulated ionization fields. We incorporate here the ef- 
fect of redshift space distortions, taking into account the 
simulated peculiar velocity field. On large scales, lin- 
ear infall boosts the spherically averaged 21 cm redshift 
power spectrum relative to its real space analogue, analo- 
gous to the 'Kaiser effect' in galaxy surveys flKaiser|1987[ 
Bharadw aTfc Ali|2004| |McQuinn et al.|20 05 ; Barkana & 
Loeb||2005p . By spherically averaging the signal we lose 
information about the ionizing sources as we ll as cosmo- 
logies , ! para meters, as was discussed e.g. in |Barkana fc 



; 65.6 Mpc/h. 



Loeb (2005). However, the first generation 21cm exper- 
iments will be sensitive mainly in the frequency direc- 
tion and have difficult y measuring the full a ngular de- 
pendence of the signal ( McQuinn et al.|[2005 ). 

The result of our power spectrum calculation is shown 
in Figure [9] for three different redshifts during reioniza- 
tion. The results are qualitatively similar to those of 
Figure[5j and can be roughly understood by decomposing 
the 21 cm power spectrum into three constituent pieces 
(FZH04): 

A» x (fc) =T 2 b [AUk) ~ lx H A 2 xS (k) + ™x 2 H A 2 ss (k)} . 

(3) 

10 We note the possibility that feedback effects, which have note 
been included in our simulations, might suppress the formation of 
the lowest mass sources an d lead in extreme c ases to a morphology 
that resembles that seen in llliev et al.l 420051) . 
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Fig. 8. — Dependence of reionization morphology on source density. In the left panel we show the ionization field from our halo-smoothing 
procedure using only sources (white points) with mass larger than M > 4 X 10 Mq (note that some sources contributing to the ionized 
regions lie in front or be hind the thin slice shown). With this choice, the number density of sources roughly matches that of M > 2x lO 9 Af0 
sources at z ~ 14 (as in |Iliev et al.|2005| |. The center panel shows the result with our usual source prescription, indicating a significantly 
more complex morphology. Finally the right panel shows, for comparison, the analytic model with M m i n = 1O 8 M0. Each panel is at 
z = 7.26, and in each case the source efficiencies are adjusted to match x lv = 0.33. 



Here A^ x refers to the ionization power spectrum, A x5 
refers to the ionization-density cross power spectrum, 
and A^ refers to the density power spectrum. Note that, 
for illustrative purposes we ignore higher or der terms 
( |McQuinn et al.)2005[ |Furlanetto et al J2006a|) although 
their effects are included in our calculations. The numer- 
ical coefficients in this decomposition come from angle- 
averaging the redshift space power spectrum. On scales 
much larger than the size of the ionized bubbles, each 
term in this decomposition is directly proportional to 
the density power spectrum, and so the 21 cm power 
spectrum is directly proportional to the density power 
spectrum. On the other hand, on very small scales one 
would expect that the 21 cm power spectrum approaches 
the density power spectrum multiplied by the neutral 
fraction squared (and a constant factor ~ 1.87 for the 
spherically averaged redshift space case). The latter is 
shown in the thin dashed curves in the Figure. The dis- 
crepancy seen is due to the significance of higher order 
terms that were neglected in Equation |3j that in rea lity 
amount to corrections of order one ( Lidzet. aT||2006[ ) . 

These qualitative trends can be seen in Figure [5] For 
further illustration, we extrapolate our predictions to 
large scales using an analytic model hybrid simulation 
(green long-dashed lines) which we based on a Gaussian 
random field with sidelength 300 Mpc/h. 

At high redshift, where the ionized regions are small, 
the 21 cm power spectrum has the shape of the den- 
sity power spectrum. At lower redshifts, it begins to 
flatten on large scales owing to the presence of ionized 
regions, before following the shape of the density power 
spectrum again on small scales. This flattening moves to 
progressively larger scales as reionization proceeds, and 
the bubbles grow larger. Our first observational han- 
dle on the characteristic sizes of HII regions at different 
stages of reionization will likely come from measuring the 
21 cm power spectrum, and observing this flattening. In 
other work, we will explore the extent to which the size 
distribution of HII regions can be extracted from future 



radiative transfer 21 cm power spectra is even better 
than the agreement between the ionization power spec- 
tra. While the ionization field in the radiative transfer 
simulation has more small scale power than the analytic 
model ionization field, the different approaches show sim- 
ilar amounts of small scale 21 cm power. This owes to 
the small-scale dominance of the A'j s (k) term in the 21 
cm power spectrum, which overwhelms the difference in 
small scale ionization power (see Figure pi). The 21 cm 
power spectrum in each analytic scheme seems to pro- 
vide a very good approximation to the results of our full 
radiative transfer simulations. Some of the difference on 
large scales may be attributable to our limited simulation 
volume, and a convergence test with increasing boxsize 
would be informative, but we leave this to future work. 

7. CONCLUSIONS 

We have presented results from a large volume radia- 
tive transfer simulation and fast numerical schemes based 
on analytic considerations, and given a detailed com- 
parison. Our basic conclusion is that the approximate 
schemes agree remarkably well with the radiative trans- 
fer simulation. 

Future work should investigate the effect of recombi- 
nations which, we anticipate, will lea d to two primary 
modifications ( Furlanetto fc Oh||2005[ ). First, recombi- 
nations will slow down reionization by requiring more 
ionizing photons to achieve a given ionization fraction. 
This should mainly act to modify the redshift evolution 
of the ionization fraction, and not the size distribution 
of HII regions at a given ionization fraction, our main 
focus in the present work. Second, ionization fronts may 
be halted upon impacting dense clumps, where the re 



combina tion rate is very high (e.g. 

2000] iShapiro et al.||2004| [Furlanetto & Oh|2005 



Miralda-Escude et al. 
This" 



meas urements of the 21 cm power spectrum (Zahn et al. 
20061. 



T^tice that the agreement between the analytic and 



latter effect might, indeed, modify the size aistribution 
of HII regions at a given ionization fraction. However, as 
long as mini-halos are destroyed by pre -heating prior to 
reionization (e.g. Oh & Haiman 20031, estimates show 
this effect is importan t only at the tail end o f reioniza- 
tion, when x- h y > 0.77 ( |Furlanetto fc Oh|2005| , which we 
do not presently simulate! 
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Fig. 9. — The 21 cm brightness temperature power spectra in 
redshift space. The solid red, short-dashed blue and dotted pur- 
ple lines show the radiative transfer, analytic, and halo-smoothing 
power spectra, respectively. The green long-dashed lines show ex- 
trapolations of the analytic predictions to large scales. Some of the 
differences in the predictions on large scales may be attributable 
to our limited simulation volume. The redshift space 21 cm power 
spectrum approaches Pgx^jl.87 (shown in the thin dashed curve) 
on small scales. The differences seen are due to the relevance of 
higher order contributions to the 21 cm power spectrum (see up- 
coming work). 



In the future we will address these issues explicitly, 
along with other refinements to our radiative transfer 
simulations. We intend to consider a more sophisticated 



prescription for the ionizing sour ces (|Springel fc Hern 



quist 



2002 Sokasian et al. 



20031, and extend the mass 



range of our sources down to the cooling mass. It will 
be interesting to examine how sensitive the 21 cm pre- 
dictions are to the assumed pro perties of the ionizing 
sources 



Furlanetto et al. 2006a I. In particular, in 



we found that the morphology and size distribution 6T 
HII regions differs dramatically from our fiducial model 
when extremely rare, bright sources dominate. This war- 
rants further quantitative investigation. Finally, we in- 
tend to examine the effect of feedback on r eionization, 
incorporating Jeans mass suppression (e.g. Barkana fc 
Loeb 20001 iBabich fc Loeb||2005 | i Kramer et aT | | 2006[) in 



reionized regions of the 1GM (Mcyuinn et al 



wmf 



In spite of these refinements, we contend that the 
agreement demonstrated in this paper illustrates that the 
analytic models are on the right track, and provide a use- 
ful complementary tool to radiative transfer simulations. 
The approximate schemes described here are very fast, 
allowing quick coverage of a large parameter space, con- 
venient f or forecasting con straints from upcoming 21 cm 



surveys (Zahn et al. 2006). Even full radiative transfer 



simulations currently have a large number of free param- 
eters related to the efficiency of the ionizing sources, the 
escape fraction of ionizing photons, and sub-grid clump- 
ing. Our numerical schemes allow one to gauge how 
the expected signal depends on these numerous, uncon- 
strained parameters. It can also be used to investigate 
non-Gaussianities in the 21 cm signal, as advocated by 
|Furlanetto et al. (2004b I, and to construct mock 21 cm 
survey volumes, providing a useful test of data analy- 
sis procedures, which are presently still under develop- 
ment. This is particularly relevant given that surveys 
like the MWA will be done in large volumes of several 
co-moving cubic Gigaparsecs, prohibitive for current ra- 
diative transfer simulations, but manageable with ana- 
lytic calculations. Finally, it might be interesting to cou- 
ple the fast analytic model schemes with a gas-dynamical 
calculation to investigate the impact of reionization on 
galaxy formation. 
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APPENDIX 
PHOTON CONSERVATION IN OUR APPROXIMATE SIMULATION SCHEMES 

The objective of this Appendix is to show that the pure FZH04 model conserves photons, but that our numerical 
schemes do not precisely conserve photons. We then discuss the implications of this finding. In the pure FZH04 
model, we can prove that the global ionization fraction is given by x — C X /coll.- This is just a reflection of photon 
conservation: as we sum up the total ionized mass from individual HII regions, no photons are lost or gained in our 
accounting of the net ionized mass. 

A rigorous proof proceeds as follows. For simplicity, we outline this proof using the pure FZH04 barrier, but the 
proof can be easily generalized to the barrier of Equation pi. Let us consider random walks in the (<5, a 2 ) plane (e.g. 
Bond et al. 1991), generated using top-hat smoothing in ftnspace. We consider the first up-crossing distributions for 
two types of barriers. First, we examine the probability that a random walk crosses the 'bubble barrier', representing 
the critical density threshold for a region to self-ionize (see Figure 1 of FZH04) . We denote the differential probability 
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Fig. 10. — 21 cm brightness temperature fluctuations. We compare 21 cm maps from the radiative transfer simulation and numerical 
scheme at three different redshifts. Each map is 65.6 Mpc/h on a side, and 0.25 Mpc//i deep, comparable to the frequency resolution of 
planned experiments, and shows a different cut then Figure [3| The ionized fractions are x\ y = 0.11,0.33 and 0.52 for z = 8.16,7.26 and 
6.89 respectively. Left column: Radiative transfer calculation with ionizing sources (blue dots). Middle column: Halo-smoothing procedure 
(see text) with sources/halos from the N-body simulation. Right column: Constant mass-to-light ratio version of FZH04, based purely on 
the initial, linear dark matter overdensity. 
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that a random walk crosses this barrier, at a resolution between a 2 and a 2 + da 2 , by dPb/da 2 . Next, we consider 
the ordinary Press-Schechter barrier, representing the critical overdensity for a region to collapse and form a halo. 
The differential probability distribution for a random walk to cross the 'collapse barrier', at a resolution between a' 2 
and a' 2 + da' 2 , is denoted by dP c /da' 2 . Similarly, the probability distribution for collapse in a region with large-scale 
overdensity 5b, on smoothing scale a 2 , is denoted by dP c (a' 2 \Sb, a 2 ) /da' 2 . The total ionized mass in a region of large 
scale overdensity 5b, at a smoothing scale a 2 , is then given by 

dMUM h ^ dP ^yi (Al) 

S dM h da 12 K ' 

Note that the conditional probability distribution in this formula is calculated by considering the fraction of random 
walks, originating at (Sb, a 2 ), t hat cross the collapse barrier at higher resolution (Lacey & Cole 1993). The mass 



calculated using Equation (Al) is precisely the ionized mass in an HII region that crosses the 'bubble barrier' at 
the point (5b, a 2 ). In order to find the tot al io nized mass in all HII regions, we merely need to integrate over all 



such crossings, i.e., we integrate Equation (Al) over a 2 weighted by the probability of crossing the bubble barrier. 
Symbolically, the total ionized mass in the IGM is then given by 

dP ^ 2 ) fji, „, da' 2 dP c (a' 2 \5 b ,a 2 ) 



d a 1 — ^V^ / dM h C M h ^ ^ — '- . (A2) 

da 2 J S dM h da 12 V ' 

This is one expression for the total ionized mass in the IGM, obtained by summing the ionized mass in all individual 
HII regions. Our proof of photon conservation is completed by showing that this 'local' expression matches a separate 
expression, proportional to the global collapse fraction. The total mass in halos is simply 

and the total, photon-conserving, ionized mass is just C, times this expression. Now, this expression, proportional to 
the global collapse fraction follows by considering the crossing distribution of the collapse barrier, i rresp ective of when 



each random walk crosses the bubble barrier. This result clearly must match that of Equation (Al) since for two 
random variables, x and y with probability distributions P(x) and P(y), J dyP(y) J dxxP(x\y) = J dxxP(x), i.e. in 
one case we are integrating ('marginalizing') over 'bubble crossings', and in the other case we are not. This proves 
that the pure FZH04 model conserves photons, and our numerical implementation of the FZH04 model with a sharp 
fc-space filter indeed conserves photons. 

In practice, however the hybrid scheme of 33] smoothes the density field with a top-hat in real space, rather than a 
sharp fc-space filter. In this case photon conservation is not guaranteed. Specifically, the expression in Equation (fl]) 
of q3] is rigorously equal to the collapse fraction only for sharp fc-space filtering, and not for real-space smoothing (see 
also McQ uinn et al.||2005[ ). One option would be to simply apply our algorithm with a sharp fc-space filter, but we 
find that this produces artificial features in our ionization maps (ringing in configuration space) . For this reason, we 
prefer to apply our algorithm using a top-hat in real space. In practice this leads to photon non-conservation at the 
20% level, with our algorithm systematically under-shooting the expected ionization, x — C, X / co n.. To compare with 
the radiative transfer simulation, we simply boost the ionizing efficiency to make up for this photon loss, matching the 
(volume-weighted) ionization fraction in the radiative transfer simulation. 

Is photon-conservation fulfilled in our improved 'halo-smoothing' scheme? We consider a simple toy problem to 
illustrate that our improved scheme also does not quite conserve photons. Imagine two equal luminosity sources in 
a uniform density field. When the ionized regions surrounding these sources begin to overlap, the spherical top-hat 
criterion can lead to somewhat unphysical features. This is sketched in the left panel of Figure TTJ Our algorithm 
does not allow for flux from one source to expand the HII region surrounding the second source. Instead of both HII 
spheres (with initial radius r-y) growing further dur ing overlap, a new ionized region arises between them, the overlap 



of two spheres with radius r 2 = 2 1 ' 3 r 1 . In Figure 11 we plot the ratio of the ionized volume in our scheme, to the 
expected, photon-conserving ionized volume. The figure clearly illustrates that our scheme generally loses photons 
as two bubbles 'merge'. The precise level of photon loss in our 'halo-smoothing' scheme will depend on the ionized 
fraction, the size distribution of the HII regions, the luminosity and bias of the sources interior to merging bubbles, 
and the rate of merging bubbles. In practice, the level of photon non-conservation in our halo-smoothing scheme is 
also at the 20% level. Again our solution is to uniformly boost the ionizing efficiency of our sources to match the 
(volume-weighted) ionization fraction in the radiative transfer simulation. Ideally, we would only boost the efficiency 
in recently merged bubbles where we expect photon loss. In practice, any error associated with this approximation 
appears small, although the higher frequency of bubble mergers at late stages of reionization makes our scheme slightly 
less reliable in this regime (see Figure ro|. 
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